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tion  for  FREDDEX  was  minimal.  A  typical  value  for  the  standard  deviation  of  the  array  was 
0.2  m.  The  lowest  spatial  mode  (the  bow  shape)  piedominated  and  cycled  in  about  16  min. 
Acoustic  degradation  was  evaluated  for  several  frequencies  and  was  always  well  below  one  dB. 
Increased  sidelobes  were  apparent,  but  they  were  still  35  dB  below  the  maximum  response. 
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A  CONSTRUCTIVE  ALGORITHM  FOR 
DETERMINING  ARRAY  SHAPE  WITH  APPLICATION 

TO  FREDDEX 

INTRODUCTION 

Background  on  the  Array  Shape  Problem 

Long,  towed  hydrophone  arrays  are  useful  in  a  sonar  context  as  high-gain,  high-resolution 
receivers.  However,  conventional  beamforming  assumes  the  array  is  straight  and  horizontal,  which  is 
generally  not  the  case  for  the  flexible,  oil-filled  hoses  presently  envisioned  for  surveillance  applications. 
Future  work  will  emphasize  computer  compensation  of  the  array  shapes  prior  to  final  beamforming. 
Recent  work  in  this  area  has  concentrated  on  measurement  techniques,  full-scale  and  model  experi¬ 
mentation,  and  analysis  of  shapes.  A  general  conclusion  is  that  only  a  few  of  the  lowest  order  spatial 
modes  have  a  significant  effect  on  array  performance  between  30  and  300  Hz. 

An  excellent  summary  of  measurement  techniques  was  presented  by  Martin  ll].  He  covered  pas¬ 
sive  and  active  acoustic  methods  as  well  as  non-acoustic  sensor  methods.  He  mentioned  analysis  tech¬ 
niques  and  their  associated  claims  of  statistical  significance.  In  general,  most  of  the  combined 
measurement/ analysis  procedures  claim  a  shape  accuracy  of  0.1  to  0.5  m.  The  actual  results  on  a 
variety  of  arrays  are  contained  in  the  classified  literature. 

Background  on  FREDDEX 

FREDDEX,  an  integrated  environmental-acoustic  experiment  conducted  in  the  Atlantic  Ocean 
near  Bermuda  in  the  spring  of  1979,  involved  several  ships,  aircraft,  and  the  use  of  satellites  on  a 
variety  of  underwater  acoustic  and  environmental  missions.  The  main  objectives  were  to  locate  and 
describe  a  large  oceanographic  thermal  anomaly  and  to  conduct  a  series  of  acoustic  measurements.  The 
underlying  purpose  of  the  experiment  was  to  assess  the  impact  of  ocean  eddies  on  the  acoustic  perfor¬ 
mance  of  present  and  future  undersea  surveillance  systems. 

In  recent  years  features  with  strong  thermal  gradients,  such  as  the  Gulf  Stream  and  warm  eddies, 
have  been  reported  as  affecting  subsurface  acoustic  signal  transmission.  These  features  have  been 
shown  by  experiment  to  reduce  signal  levels  by  as  much  as  15  dB,  depending  on  geometry,  source  and 
receiver  depths,  and  "strength"  of  thermal  gradients.  FREDDEX  was  designed  to  extend  previous 
knowledge  to  a  cold  eddy  environment  and  to  expand  the  analysis  to  include  high-resolution  measure¬ 
ments  of  the  performance  of  towed  arrays.  Measurements  such  as  array  signal  gain,  resolving  power, 
and  target  bearing  accuracy  were  included.  The  initial  acoustic  results  from  a  time-limited  data  set  indi¬ 
cate  good  towed-array  performance  for  the  particular  oceanographic  feature  studied  and  the  specific 
array  used. 

To  conduct  a  thorough,  broadly  supportive  exercise,  ten  different  organizations  collaborated  to 
provide  expertise  in  the  various  source,  receiver,  and  environmental  areas  of  concern.  Some  experi¬ 
ments  conducted  in  support  of  FREDDEX  were:  satellite  and  over-the-horizon  radar  remote  sensing  of 
ocean  surface  characteristics,  airborne  sensing  of  shallow  thermal  structure,  shipboard  in-situ  measure¬ 
ments  of  deep  thermal  structure,  and  towed  and  bottom-mounted  receiver  measurements  of  sound 
transmitted  from  moored  and  towed  acoustic  projectors. 
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Description  of  the  FREDDEX  Array  Shape  Experiment 

One  particular  receiver  system  used  during  the  three  weeks  at  sea  phase  of  FREDDEX  was  a 
long,  flexible  hydrophone  array  towed  by  a  surface  ship  at  depths  of  several  hundred  meters.  The  array 
consisted  of  60  hydrophones  spaced  10  m  apart.  This  corresponds  to  32  wavelengths  at  83  Hz  or  88 
wavelengths  at  227  Hz.  The  array  performance  depends  on  many  factors  including  the  spatial  coher¬ 
ence  of  the  incon.ing  acoustic  wavefronts  and  the  towed  receiver  shape  and  dynamics.  The  various  fac¬ 
tors  influencing  performance  confound  the  results  and  thus  are  examined  independently  ;o  order  to 
determine  their  relative  contributions, 

A  six  hour  time  period  at  sea  was  devoted  exclusively  to  measurements  of  array  shape  and  array 
dynamics.  The  array  deformation  experiment  consisted  of  the  receiver  ship  towing  the  hydrophone 
array  at  a  depth  of  180  m  for  three  hours  and  260  m  for  an  additional  three  hours  at  a  position  essen¬ 
tially  broadside  to  the  array.  For  these  six  hours  the  acoustic  source  towing  ship  moved  along  a  track 
parallel  to  the  receiver  towing  ship  at  a  range  of  about  2.8  km.  The  acoustic  projector  was  towed  at  a 
depth  of  200  m.  A  7  ms  pulse  (2  cycles  at  300  Hz)  with  ISO  Hz  bandwidth  was  transmitted  from  the 
acoustic  projector  every  15  s.  The  puise  received  by  the  array  was  on  the  order  of  3  cycles  due  to 
source/ transducer  ringing.  Absolute  arrival  times  to  each  hydrophone  were  obtainable  by  providing 
each  ship  with  a  crystal  controlled  clock  (synchronized  within  1  ms  to  a  common  radio  source, 
(WWV))  thus  allowing  accurate  determination  of  the  lowest  order  spatial  mode,  i.e.  a  bow  shape,  as 
well  as  the  finer  scale  shape. 

This  report  documents  an  algorithm  developed  for  determining  the  shape  of  the  towed  array  as 
well  as  the  dau  processing  techniques  used  on  the  collected  acoustic  FREDDEX  data  set  to  determine 
the  array  shape  and  its  effect  on  acoustic  performance. 

THE  ARRAY  SHAPE  ALGORITHM 

Description  of  the  Algorithm 

At  any  fixed  time  r,  let  S  denote  the  source  and  P\.P2,  — P\  denote  the  A  hydrophone  positions 
(the  geometry  is  shown  in  Fig.  1).  Put  9,  =  k  ( P ,  SP,+\ ),  <b,  =  k  ( P,P\S )  and  </>,  =  k  (P,  P\  P,-\)- 

i 

Let  a,  =  \P]P,  I  and  b  denote  the  fixed  uniform  distance  between  array  elements.  Define  9{l>  =  £  9j. 

j- i 

Assume  that  the  A  time  deLys  A Zj .  Ar2 . Ar,v  between  the  source  and  the  A  hydrophones  are 

known.  Let  c  denote  a  representative  sound  speed  in  the  ocean  area  of  operation.  Then  the  distances 
d\,di . dN  between  the  hydrophones  and  the  source  are: 

d,  =  c  At,,  i  =  I . A. 


For  fixed  /',  let  X  denote  the  intersection  of  P ,  P,  with  SP>_ |  (extended  if  necessary). 


Assuming  that  <z,  .|  is  known,  an  iterative  procedure  for  determining  the  rth  hydrophone  position 
will  be  derived.  Since  a2  is  the  distance  between  consecutive  hydrophones,  this  procedure  can  then  be 
iteratively  applied  to  determine  all  hydrophone  positions.  By  the  Law  of  Cosines: 


9,  =  arc  cos 


d,2+d,lt~b2 
2  d,dl+] 


i  =  1 


Hence  9U  21  and  9,  ,  are  known.  Again  by  the  Law  ol  Losiries: 

a,2  =  d{  +  d}  -  2</,rf,  cos  9U  " 
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Additional  Remarks  on  the  Algorithm 

Since  it  is  desired  to  use  the  array  shape  for  additional  processing,  it  is  useful  to  associate  a  coor- 

k  _ 

dinate  system  with  the  solution.  Define  d>(A)  =  £  for  3  4  k  ^  N.  Let  P\  =  (0,0)  and  P\P:  lie 

i- 3 

along  the  positive  x  axis.  Then: 

Pt  =  (0,0), 

P  2  =  (a2,0), 

Pj  =  (a3  cos  4>{1),a}  sin 
P^  =  («4  cos  <£l4l,a4  sin  <£(4)), 

Pn  =  (ff.v  cos  <t>{S).as  sin  d>(V)).  and  5  =  ( d ,  cos  (<£,-t/r3).  sin  (03-i /»,)) 
where  i/»3  is  positive  and  <63  is  taken  with  its  sign. 

Synthetic  data  were  used  to  test  the  algorithm.  The  array  was  assumed  to  be  piecewise  linear  and 
distances  from  an  acoustic  source  to  each  hydrophone  were  computed  assuming  the  perpendicular  dis¬ 
tance  from  the  acoustic  source  to  the  array  was  1.5  miles  (as  in  the  experiment).  These  distances  were 
then  input  into  a  computer  program  designed  to  implement  the  algorithm.  The  original  linear  shape 
was  successfully  reproduced.  To  test  the  algorithm  for  statistical  robustness,  individual  hydrophones 
within  the  simulated  array  were  varied  by  amounts  up  to  several  meters.  The  shapes  were  then  recom¬ 
puted  and  it  was  found  that  the  altered  point (s)  had  almost  no  effect  on  the  computed  positions  for  the 
remaining  hydrophones.  Finally,  several  piecewise  linear  shapes  were  input  into  the  algorithm,  and  it 
was  found  that  the  algorithm  reproduced  them  exactly.  Figure  2  shows  the  results  from  a  test  case  in 
which  values  for  a  three  segment  array  producing  a  deep  bow  effect  were  computed  and  input  into  the 
algorithm.  The  shape  was  correctly  reproduced. 


0  100  200  300  400  500  600 


METERS 

Fig  2  —  A  lest  data  case 
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For  FREDDEX,  depth  sensors  were  placed  at  the  front  and  the  tail  of  the  array,  and  it  was  found 
that  depth  variations  were  insignificant,  i.e.  the  vertical  tilt  was  less  than  0.5°.  Accordingly,  in  develop¬ 
ing  an  algorithm  to  process  FREDDEX  data,  depth  could  be  ignored.  In  addition  it  was  assumed  that 
the  array  spacing  was  known  and  fixed  and  that  the  array  was  piecewise  linear. 

DATA  PROCESSING 


Digitizing  and  Time  Delays 

The  next  several  sections  discuss  the  implementation  of  the  algorithm  to  the  FREDDEX  acoustic 
data.  Figure  3  traces  the  processing  and  display  steps  from  the  analog  acoustic  tapes  through  the  com¬ 
putation  of  array  shape  coordinates  for  the  experiment.  Figure  4  shows  the  flow  of  computations  made 
from  this  file  of  array  shape  coordinates.  In  both  cases,  the  program  names  are  given  along  the  solid 
lines.  FORTRAN  source  listings  for  the  major  data  processing  programs  are  given  in  the  appendix  of 
this  report. 


Fig.  3  _  Flow  chart  of  data  processing  to  obtain  array  shape  positions 


Fig.  4  —  Flow  chart  of  processing  using  file  of  array  shape  positions 
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Existing  computer  programs  were  used  for  the  digitization  and  real-time  display  of  collected  ana¬ 
log  acoustic  data.  Twenty-seven  alternate  channels  of  data  spaced  at  20  m  were  digitized  at  a  sampling 
rate  of  4096  samples  per  second.  This  produced  eleven  samples  per  cycle  at  the  highest  frequency  of 
the  pulse  (375  Hz).  Three  digital  data  tapes  were  produced  from  the  6  hours  of  data  for  further  pro¬ 
cessing.  Figure  5  is  an  example  of  the  real-time  display  used  to  monitor  the  digitizing  process.  Hydro¬ 
phones  2,  15,  and  27  are  monitored  and,  from  the  time  of  arrival,  one  can  see  that  the  closest  hydro¬ 
phone  to  the  source  was  between  hydrophones  2  and  15  for  this  time  slice.  A  digitizing  window  length 
of  125  ms  per  channel  was  chosen  to  capture  the  largest  expected  bow  shape  of  the  array.  Hydrophone 
15  was  placed  between  31.25  and  62.5  ms  since  the  experiment  attempted  to  align  the  two  ships  so  that 
hydrophone  15  would  be  closest  to  the  acoustic  source. 


_ I _ I 


0  31  25  62  5  93  75  125 

TIME  l^sec) 


Fig .5  —  Real-time  digitization  display 


The  next  step  was  to  produce  time  delays  using  program  ARVTM  ARVTM  read  the  digitized 
tape  and  computed  time  delays  which  were  stored  on  a  disc  file  for  further  processing.  In  order  to 
eliminate  transients  which  occurred  during  the  digitizing  process,  any  detection  which  occurred  outside 
a  two  second  window  for  the  shot  was  not  processed.  For  valid  times  the  largest  negative  peak  for  each 
phone  was  found  and  a  threshold  was  set  at  1/2  the  absolute  value  of  this  peak  Next,  the  index  of  the 
first  negative  value  which  (negatively)  exceeded  this  threshold  was  determined,  and  the  first  positive 
value  of  greater  index  was  located.  The  ping  arrival  time  was  then  computed  from  the  zero  crossing  of 
a  linear  interpolation  between  the  two  points.  Once  this  had  been  done  for  all  27  hydrophones,  the 
sequence  of  indicies  corresponding  to  each  negative  peak  position  was  examined  to  assure  that  it  had 
only  one  turning  point.  If  an  erroneous  peak  was  located,  the  average  of  the  indices  of  the  peaks  of  the 
two  surrounding  hydrophones  was  used  as  an  estimate  of  the  peak  and  the  zero  crossing  was  recom¬ 
puted.  Finally,  time  delays  were  obtained  for  each  hydrophone  by  using  1  >  the  header  record  to  obtain 
the  time  delay  to  the  start  of  the  digitized  data,  2)  the  index  of  the  Iasi  negative  point  before  the  zero- 
crossing  to  obtain  the  time  to  this  point  from  the  start  of  the  data,  and  3)  a  linear  interpolation  between 
the  last  negative  and  first  positive  value  to  obtain  the  time  of  the  zero  crossing  This  detection  process 
proved  quite  successful  allowing  only  a  few  instances  of  detection  errors,  almost  all  of  which  occurcd  in 
the  first  90  minutes  of  data  prior  to  the  source  level  being  increased 
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Table  !  —  Typical  Time  Delay  Calculation 


START 

Tirt: 

171  6 

6 

1 

PHQNt 

MAX 

NEGPK 

2ER0+ 

+  VAL 

-  VAL 

TT1  DELAY 

1 

8064 

ill 

2i7 

1406 

>816 

1.834074 

2 

9856 

207 

212 

3104 

-112 

1.032772 

3 

9920 

203 

208 

1568 

-1216 

1.831894 

4 

8544 

199 

203 

944 

-2256 

1.830739 

5 

9008 

196 

201 

128 

-2752 

1.030311 

6 

10352 

193 

198 

2528 

-224 

1.829366 

7 

9600 

(90 

196 

2096 

-112 

1.828870 

8 

9056 

188 

193 

512 

-2368 

1.828326 

9 

10080 

187 

191 

480 

-2480 

1 . 827841 

10 

7952 

185 

190 

1040 

-1456 

1.827535 

1 1 

10096 

184 

189 

1904 

-1568 

1.027259 

12 

8608 

184 

188 

160 

-3056 

1.827136 

13 

94 40 

183 

188 

640 

-2432 

1.827098 

14 

10432 

183 

188 

2368 

-1152 

1.826984 

15 

8992 

185 

190 

1552 

-99? 

1 . 0274AA 

16 

9646 

185 

191 

2432 

-256 

1.027660 

17 

7632 

188 

193 

1792 

-304 

1.828168 

10 

10272 

190 

194 

464 

-3000 

1.828581 

19 

7152 

191 

197 

1792 

-80 

1.029112 

20 

9952 

195 

201 

1920 

-720 

1.830145 

21 

9344 

198 

203 

560 

-2480 

1.830766 

22 

6592 

202 

209 

1360 

-128 

1.632052 

23 

6480 

206 

211 

1488 

-1152 

1.832626 

24 

10592 

211 

216 

168 

-2592 

1.033970 

25 

9488 

216 

221 

1312 

-1520 

1.835092 

26 

7800 

223 

228 

1 488 

-032 

1.836757 

27 

8912 

228 

233 

64 

-2496 

1.838129 

Table  1  illustrates  the  process  of  producing  time  delays  for  a  fixed  time  slice  (time  =  171  6  6). 
For  each  hydrophone  the  largest  negative  peak  was  found.  The  first  negative  peak  exceeding  ihe  thres¬ 
hold  for  each  hydrophone  was  then  computed  and  listed  in  the  column  labeled  NEGPK  in  terms  of  its 
index  within  the  time  segment,  r  e.  from  I  to  512.  The  column  labeled  ZUR 04-  contains  the  index  of 
the  first  positive  value  following  the  first  negative  peak  and  +V,\L  contains  that  value.  The  column 
labeled  -VAL  contains  the  value  of  the  previous  index.  Interpolation  between  +VAL  and  -VAL  then 
allowed  accurate  computation  of  the  time  delay  as  seen  in  column  labeled  TM  DELAY.  The  time 
delays  are  then  stored  in  a  data  file  on  disc  for  further  processing. 


Program  MRGTM  then  merged  the  three  time  delay  files  into  a  single  file  for  further  processing. 
At  the  same  time  a  further  check  was  made  to  assure  that,  for  a  given  time,  the  ping  appeared  within 
the  proper  time  interval.  Using  estimated  arrival  times  derived  from  the  range  and  sound  speed  pings 

which  were  not  detected  between  s  and  s  after  the  quarter  minute  marks  were  eliminated. 


Array  Shapes 

Next,  time  delays  were  used  to  produce  array  shapes  using  program  ARSHP,  which  implements 
the  algorithm  previously  described  and  produces  plots  of  single  scan  array  shapes  as  well  as  water  fall 
plots  of  array  shapes  over  user  defined  time  periods.  Using  the  measured  sound  speed  of  1522  m/s, 
time  delays  were  converted  into  distances  and  the  array  shape  algorithm  was  applied,  resulting  in  a  set 
of  x-y  coordinates  for  each  time  slice.  In  order  that  successive  shapes  could  easily  be  prepared  and 
beamforming  could  be  performed,  a  best-fit  least-squares  line  was  calculated  for  the  shape  at  each  fixed 
lime  and  the  x-y  coordinates  were  translated  and  rotated  by 

x,  =  x,  cos  (a)  +  O',  —  /3)  sin  (a), 
y,  =  -x,  sin  («)  +  O,  -  p)  cos  («), 

where  a  is  the  angle  that  the  least  squares  fit  line  makes  with  the  x-axis  and  (3  is  the  intercept  with  the 
y-axis.  Plotting  routines  were  then  used  to  get  individual  snapshots  of  the  array  and  to  produce  water¬ 
fall  plots  of  array  shape  as  a  function  of  time.  Figure  6  is  an  example  of  a  single  snapshot  of  the  array. 
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Eig  6  —  Typical  array  sh».pe 


This  snapshot  was  taken  about  five  hours  into  the  experiment  (at  time:  171  6  0)  and  illustrates  the 
hydrophone-to-hydrophone  variability  found  by  applying  the  algorithm  without  averaging.  Larger  varia¬ 
tions,  as  seen  at  hydrophone  22,  have  often  been  found  to  be  consistant  over  a  small  time  interval,  i.e. 
a  variation  of  about  0.4  m  from  an  expected  position  would  occur  for  a  15-  to  30-minute  time  period 
and  the  hydrophone  would  fall  back  into  position.  This  deviation  corresponds  to  an  error  of  approxi¬ 
mately  one  sample  and  is  too  small  to  be  a  peak  detection  problem.  One  possible  explanation  would  be 
a  jump  of  one  position  in  the  digitizer. 

The  next  step  in  the  processing  was  to  correct  for  hydrophone-to-hydrophone  phase  delays.  To  do 
this,  the  section  of  data  between  5:00  and  6:30  was  selected  because  the  averaged  shape  for  this  time 
period  exhibited  only  linear  characteristics.  Time  delays  which  would  linearize  the  shape  were  com¬ 
puted  for  each  hydrophone.  Table  2  shows  the  absolute  time  delays  and  Table  3  shows  the  phase 
differences  (relative  to  hydrophone  14).  The  average  array  for  the  entire  seven  hours  of  data  is  shown 
in  Tig  7  The  smaller  deviations  noted  around  the  1/4  and  3/4  points  on  the  array  are  accounted  for 
by  the  predominance  of  the  first  order  mode  in  the  array.  The  slight  bow  in  the  array  is  due  primarily 
to  the  initial  12  minutes  of  data  when  the  ships  were  still  maneuvering  to  correctly  align  themselves  for 
the  experiment  Next,  program  ARSH7  was  used  to  process  the  time  delay  data  file  as  previously  dis¬ 
cussed  with  the  phone-to-phone  phase  corrections  and  to  store  the  results  in  a  new  data  file  containing 
.v-.v  coordinates  of  the  array  shapes  over  the  entire  data  set.  The  waterfall  plot  program  was  then 
modified  to  work  on  the  x-y  coordinates.  Figure  8  shows  a  waterfall  plot  of  array  shape.  The  shapes 
have  been  averaged  in  two  minute  intervals  and  the  display  contains  one  hour  of  data.  The  sinusoidal 
oscillation  of  the  array  as  a  function  of  time  is  easily  seen  with  a  period  of  about  16  minutes. 

A  very  small  number  of  bad  data  points  still  existed  in  the  data  set  at  this  point  for  a  variety  of 
reasons.  For  example,  when  the  sequence  of  indexes  of  negative  peaks  was  searched  for  nonmonotonic 
values  corresponding  to  false  peak  detection,  no  check  was  made  to  see  if  the  minimum  value  was 
correctly  detected.  To  eliminate  these  problems,  hydrophones  which  differed  by  more  than  one  meter 
from  both  surrounding  hydrophones  were  replaced  by  the  average  y-coordinate  value  of  the  two  sur¬ 
rounding  hydrophones.  This  resulted  in  editing  about  60  points  out  of  a  total  of  over  32,000.  In  addi¬ 
tion,  it  was  noted  that  on  hydrophone  1,  the  same  values  would  reoccur  on  occasion,  once  for  23  con¬ 
secutive  times.  This  was  clearly  due  to  some  error  in  the  digitization,  most  likely  a  problem  with  read¬ 
ing  the  first  channel  on  the  tape.  In  those  instances  where  the  values  were  identical  for  successive 
times,  the  first  time  delay  was  replaced  by  the  extrapolated  value  obtained  from  hydrophone  2  and 
hydrophone  3.  A  check  was  made  of  all  channels  to  see  if  this  problem  had  occurred  elsewhere,  but  it 
had  not. 
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Table  2  —  Phase  Delays  Relative  to  Hydrophone  14 


FREQ 

PHONE 

48.5 

r 

83.3 

— 

111.6 

— 

162.0 

- 

224.0 

227.0 

230.0 

235.0 

350.0 

1 

-l.i 

-1.8 

-2.5 

-3.6 

-4.9 

-5.0 

-5.1 

-5.2 

-7.7 

2 

-2.6 

-4.4 

-5.9 

-8.6 

-11.9 

-12.1 

-12.3 

-12.5 

-18.6 

3 

-2.7 

-4.7 

-6.3 

-9.1 

-12.6 

-12.7 

-12.9 

-13.2 

-19.7 

4 

-1.9 

-3.2 

-4.3 

-6.2 

-8.6 

-8.7 

-8.9 

-9.1 

-13.5 

5 

-5.0 

-8.6 

-11.5 

-16.7 

-23.1 

-23.5 

—  23.8 

-24.3 

-36.2 

6 

-0  9 

-1.6 

-2.1 

-3.1 

-4.3 

-4.3 

-4.4 

-4.5 

-6.7 

7 

-3.0 

-5.2 

-6.9 

-10.0 

-13.9 

-14.1 

-14.2 

-14.6 

-21.7 

8 

-2.1 

-3.6 

-4.9 

-7.1 

-9.8 

-9.9 

-10.0 

-10.2 

-15.2 

9 

-2.5 

-4.3 

-5.7 

-8.3 

-11.5 

-11.7 

-11.8 

-12.1 

-10.0 

10 

-4.3 

-7.4 

-14.5 

-20.0 

-20.3 

-20.5 

-21.0 

-31.2 

11 

-3.9 

-6.6 

-8  9 

-12.9 

-17.8 

-10.1 

-10.3 

-10.7 

-27.0 

12 

-3.2 

-5.5 

-7.4 

-10.7 

-14.8 

-15.0 

-15.2 

-15.5 

—23.1 

13 

-2.7 

-4.6 

-8.9 

-12.3 

-12.5 

-12.7 

-12.9 

-19.3 

14 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

15 

-2.3 

-3.9 

-5.2 

-7.6 

-10.5 

-10.6 

-10.8 

-11.0 

—  16.4 

16 

-0.8 

-1.4 

-1.9 

-2.7 

-3.8 

-3.8 

-3.9 

-4.0 

-5.9 

17 

-1.8 

-3.1 

-4.1 

-5.9 

-8.2 

-8.3 

-8.4 

-8.6 

-12.9 

18 

-4.3 

-7.5 

-10.0 

-14.5 

—20.1 

-20.3 

-21.1 

-31.4 

19 

-0.9 

-1.5 

-2.0 

-4.1 

-4.2 

-4.2 

-4.3 

-6.4 

20 

-3.0 

-5.2 

-7.0 

-14.0 

-14.1 

-14.3 

-14.6 

-21.8 

21 

-0.3 

-0.5 

-0.7 

-1.5 

-1.5 

-1.5 

-1.5 

-2.3 

22 

-4.3 

-14.4 

-19.9 

-20.2 

-20.5 

-20.9 

-31.1 

23 

-1.8 

-6.1 

-8.5 

-8.6 

-8.7 

-8.9 

-13.2 

24 

-3.0 

-7.0 

1 

— 14.0 

-14.2 

-14.4 

-14.7 

-21.9 

25 

0.0 

0.0 

0.0 

mm 

0.1 

0.1 

0.1 

0.1 

26 

-4.8 

-8.2 

-11.0 

1E1S51 

-22.1 

-22.4 

-22.7 

-23.2 

-34.5 

27 

-3.0 

-5.1 

-6.8 

-9.9 

-13.7 

-13.9 

-14.1 

-14.4 

-21.4 

Table  3  —  Absolute  Time  Delays  (ms)  for  each  Hydrophone 


1 

2 

3 

4 

5 

6 

7 

8 

9 

0.000 

-0.007 

-0.015 

0.034 

-0.146 

0.031 

-0.031 

HiMIil 

-0.002 

10 

11 

12 

13 

14 

15 

16 

17 

18 

-0.107 

-0.080 

-0.042 

-0.012 

MBM 

0.011 

0.094 

19 

20 

21 

22 

24 

25 

26 

27 

0.090 

-0.032 

-0.123 

-0.106 

0.036 

-0.033 

0.148 

-0.133 

-0.029 
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PHONE 

Fig.  7  —  Array  average  for  the  entire  event 
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Statistics 

Now  lhal  ihe  array  positions  were  determined,  it  was  possible  to  look  at  array  statistics  and  to  do 
beamforming.  First,  throughout  the  processing  it  had  been  assumed  that  the  array  elements  were  20  m 
apart.  Program  XCHK1  examined  ihe  computed  .v  positions  (after  rotation)  and  verified  that  variations 
from  multiples  of  20  were  indeed  insignificant. 

Next,  program  ARSTA  computed  the  mean  and  standard  deviation  (sigma)  for  the  array  averaged 
over  arbitrary  (user  defined)  time  intervals.  Figures  9a  and  9b  are  examples  of  this  type  of  plot.  Six¬ 
teen  minutes  of  data  centered  about  1713  0  are  seen.  The  full  oscillation  of  the  bow  shape  is  again 
observed.  The  maximal  bow  of  approximately  I  m  seen  at  both  171  2  58  and  171  3  4  is  typical.  To 
further  examine  the  array  variability,  plots  were  made  of  sigma  as  a  function  of  time  for  user  input 
averaged  intervals.  Figure  10  shows  this  type  of  plot  for  a  two-minute  average.  Sigma  is  seen  to  run  at 
about  0.2  m  except  at  times  when  the  array  has  maximal  sinusoidal  peaks.  The  gaps  in  Figure  10  occur 
at  times  when  the  acoustic  source  was  not  operational.  Figure  1 1  shows  plots  of  sigma  bins  vs  number 
of  occurrences  for  bin  size  0.02  m.  This  plot  shows  that  0.2  m  is  the  typical  sigma  over  the  experiment. 
Since  0.2  m  is  small  compared  to  the  shortest  wavelength  in  FREDDEX  (approximately  4  m),  it  is  clear 
that  the  array  deformation  caused  by  towing  had  an  insignificant  acoustic  effect. 


to 

cc 


5 


c n 


cc 


2 


to 

(T 


Ui 

2 


C/7 

CC 


2 


HYDROPHONE 


HYDROPHONE 


Fig  9a  —  Two-minute  averages  of  array  shapes 


Fig.  9b  —  Two-minute  averages  of  array  shapes 
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Fig.  10  —  Sigma  of  (wo-minute  averages  over  the  entire  event  Fig.  1 1  —  Distribution  of  sigma  for  a  bin  size  of  0.02  m 


Beamforming 

Beam  patterns  were  produced  from  the  file  of  x-y  coordinates  by  program  BMFRM.  Beam  pat¬ 
terns,  averaged  over  any  user  defined  time  intervals,  were  computed  for  arbitrary  (user  defined)  angles 
of  arrival  (<£)  and  frequency.  The  27-hydrophone  subarray  used  to  compute  shapes  was  both  interpo¬ 
lated  and  extrapolated  to  reproduce  the  original  60  hydrophone  array  spaced  10  m  apart.  A  Hanning 
window  was  applied  and  a  complex  FFT  was  taken.  The  array  was  then  reordered  (since  the  FFT 
results  are  not  ordered  by  angle)  and  normalized  and  the  power  was  taken  to  compute  the  beam  pat¬ 
tern.  The  3-dB  down  points  and  width  were  then  computed  and  a  plot  of  computed  and  nominal  beam 
patterns  produced.  Figure  12  shows  the  average  beam  pattern  for  the  entire  test  computed  for  a  per¬ 
pendicular  arrival  path  (<£  =  0°)  and  a  frequency  of  227  Hz.  Figure  13  shows  a  typical  beampattern 
over  a  15-minute  average  (start  time:  171  3  0)  for  <t>  =  0°  and  a  frequency  of  83  Hz.  Figure  14  shows  a 
one-hour  averaged  beam  pattern  for  362  Hz  (start  time.  171  3  0).  In  these  figures,  the  actual  array 
represents  the  beam  pattern  computed  for  the  previously  obtained  array  shapes  and  then  averaged  over 
the  stated  time  interval.  The  linear  array  beam  pattern  denotes  the  beam  pattern  that  would  be 
obtained  if  there  were  no  spatial  distortion  of  the  array,  i.e.  the  array  is  linear.  As  would  be  expected 
from  the  minimal  spatial  distortion  which  the  array  was  found  to  have  during  towing,  no  significant 
degradation  of  array  performance  is  noted.  In  Fig.  12  the  array  has  a  peak  degradation  of  0.11  dB  and 
the  3-dB  down  beamwidth  is  0.92  degrees.  In  Fig.  13  the  peak  degradation  is  0.013  dB  and  the  3-dB 
down  beamwidth  is  2.51  degrees.  Figure  14  has  a  peak  degradation  of  0.29  dB  and  3-dB  beamwidth  of 
0.59  degrees.  Large  increases  in  sidelobe  levels  are  seen  in  these  figures;  however,  these  sidelobes  are 
more  than  35  dB  down  and  thus  will  have  no  effect  upon  the  acoustic  analysis  problem. 

CONCLUSIONS 

1.  The  algorithm  worked  very  successfully  on  the  FREDDEX  data  set.  Small  perturbations  in 
the  array  were  successfully  reconstructed  by  the  algorithm.  The  high-source  levels  used  during  the 
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Fig.  12  —  Averaged  beam  pattern  at  227  Hz  for  the  entire  event 
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Fig  13  —  Averaged  beam  pattern  at  83  Hz  (15-minute  average) 

experiment  were  essential  to  allow  the  accurate  computation  of  arrival  limes  required  by  the  algorithm. 
Interpolation  of  the  zero  crossing  was  required  to  compute  the  arrival  time  with  sufficient  accuracy  for 
the  algorithm.  Also  the  range  was  close  enough  and  the  source  and  receiver  depths  deep  enough  to 
resolve  direct  and  surface  reflected  paths. 

2.  For  FREDDEX,  there  was  little  array  deformation.  The  computed  typical  deviation  of  0.2  m 
is  small  compared  to  the  shortest  wavelength  in  FREDDEX  of  approximately  4  m.  These  results  are 
not  typical  since  the  array  was  being  towed  at  the  rather  high  speed  of  5  knots.  However,  the  algorithm 
would  be  expected  to  work  well  for  data  with  greater  array  deformations. 


13 


ROSENBLUM  AND  DEL  BALZO 


Fig  14  —  Averaged  beam  pattern  at  362  Hz  (one-hour  average) 


3.  The  lowest  spatial  mode  (the  bow  shape)  cycles  is  about  16  minutes.  No  significant  effect 
from  higher  order  modes  was  seen. 

4.  Since  two  ships  are  needed,  current  techniques  for  acquiring  this  type  of  data  are  very  expen¬ 
sive.  Other  potential  techniques  such  as  a  distribution  of  compasses  along  the  array  or  using  a  far-field 
beacon  source  for  bearing  and  shape  computation  merit  investigation. 
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APPENDIX 


This  appendix  contains  FORTRAN  source  listings  of  the  major  programs  used  in  the  data  process¬ 
ing.  The  reader  should  refer  to  Figs.  3  and  4  where  flow  charts  of  the  programs  and  their  functions  are 
given.  The  programs  were  written  to  operate  on  a  Digital  Equipment  Corporation  (DEC)  11/45  com¬ 
puter  system  operating  under  RSX-1 1M. 


PROGRAM  ARVTM 


THIS  PROGRAM  WILL  PROCESS  THE  DIGITIZED  DATA  TAPES  AND 
PRODUCE  A  FILE  OF  TIME  DELAYS  FOR  EACH  PHONE.  A  SECOND  FILE 
WITH  MAX  AND  MIN  FOR  EACH  PHONE  WILL  ALSO  BE  PRODUCED. 

THE  ARRIVAL  TIME  IS  DEFINED  AS  THE  1ST  ZERO  CROSSING 
AFTER  THE  1ST  LARGE  NEGATIVE  PEAK.  THIS  IS  DETERMINED 
BY  FINDING  THE  MAXIMUM  NEGATIVE  VALUE  AND  THEN  THRESH- 
HOLDING  AT  \/2  THIS  VALUE.  INTERPOLATION  IS  PERFORMED 
TO  LOCATE  THE  ZERO  CROSSING  AND  SOME  CHECKS  ARE  MADE 
TO  ALLEVIATE  PROBLEMS  DUE  TO  TRANSIENTS  AND  BAD  MEASURE¬ 
MENTS  . 

PROGRAMMER:  L.  ROSENBLUM  VERSION:  4-H-80 


RECORD  STRUCTURE  FOR  OUTPUT  FILE  TIMDL 


WORD  1 
l.JORD  2 
UORD  3 
WORD  4 

UORD  5/UORD  6 
LORDS  7-S0 


JULIAN  DAY 
HOUR 
MINUTE 
SECOND 

FRACTION  OF  SECOND 

TIME  DELAYS  FOR  PHONES  1-27 


RECORD  STRUCTURE  FOR  OUTPUT  FILE  MAXVL 
LORDS  1-6  TIME  AS  ABOVE 

UORDS  7-36  MAX  VALUES  FOR  PHONES  1-30 

UORDS  34-66  MIN  VALUES  FOR  PHONES  1-30 

DIMENSION  I8UFC512). ISTOP(3). IMAX(30). IPOSC30) , ITIMC6) 
DIMENSION  ICHKC30) . IPOS1 (30) > IZERPC30) . IZERN(30) ,TDEL(30) 
DIMENSION  IARR 1 (60) . IARR0C66). IMAX1 (30) . ITIM2(6) 

EQUIVALENCE  ( IARR1 ( 1 ) . ITIMC1 ) ) . ( IARR I C7) , TDEL) . 

I  ( IARR0< 1 ) . ITIM2 ( 1 ) ) . ( I ARR0 (7) , IMAX1 (!)).( I ARR0 (37)  ,  IMAX( 1) ) 
DATA  IC0NS.LU0UT/5.5/ 

C  GET  INPUT  PARAMETERS 

TYPE  *. '  ENTER  STOP  TIME:  HOUR, MIN. SEC' 

ACCEPT  *. ISTOP(l), IST0P(2) . IST0P(3) 

TYPE  *, '  ENTER  MAG  TAPE  UNIT  *' 

ACCEPT  *. IUNIT 
C  SET  UP  OUTPUT  FILES 

CALL  ASS IGNU.'SY: TIMDL. TP3') 
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DEFINE  FILE  1  (500. 60, U, IPNT1) 

CALL  ASS  I GN ( 3 , ' SY : MAXVL . TP3 ' ) 

DEFINE  FILE  3  (500. 66, U, IPNT0) 

IPNT1 = 1 
IPNT0= 1 

CALL  ASSIGN (2. 'SY:LR0UT'  ) 

C  READ  THE  HEADER  FILE 

10  CALL  MTREAD( IBUF.512. IER.0. IUNIT. NUMB) 

DO  5  1=1,6 
ITIM( I ) = IBUF (1+6) 

5  CONTINUE 

TYPE  *,'  TIME:  ' . ( ITIMC J) . J= 1 ,4) 

C  SKIP  TIMES  WHEN  ARRAY  UAS  OUT  OF  WATER  OR  ON  DAY  IN  OR  OUT. 

IF( IBUF(8) .EQ. I .AND. IBUF(9) .GT. 16)GD  TO  22 
IF ( IBUF (8) .EQ. 2 . AND. IBUF (9) .LT. 19) GO  TO  22 
IF( IBUF (8) .EQ.4. AND. IBUF (9) .GT.20)  GO  TO  22 
IF ( IBUF (8) . EQ .5. AND. IBUF (9) .LT. 6)  GO  TO  22 
C  CHECK  FOR  STOP  TIME 

IF ( IBUF (8) .LT. ISTOP (1)5 GO  TO  15 

IF( IBUF (8) .EQ. I STOP ( 1) .AND. I8UF(9) ,LT. I STOP (2) )GO  TO  15 
IF ( IBUF (8) . EQ . ISTOP ( 1) . AND . IBUF (9) . EQ. ISTOP (2) .AND. 

1  IBUF ( 10) .LT. ISTOP (3) ) GO  TO  15 
IF ( IBUF (7) .EQ. 171) GO  TO  99 
CALL  SKIPF( IUNIT. 1) 

GO  TO  10 

C  CHECK  THAT  PING  FALLS  IN  PROPER  SECOND  -  IF  NOT  REJECT  FILE 

15  DO  20  1  =  1.60.  15 
J  =  1+2 

IF ( IBUF (10) .GE . I .AND . IBUF ( 10) .LE. J) GO  TO  30 
20  CONTINUE 

22  UR  I TE (2,25) ( IBUF (J),J=?,10) 

25  FORMATdX. 'TIME  =  *  ,4I4.2X.'NOT  IN  BOUNDS') 

C  GO  TO  NEXT  TIME  HEADER 

CALL  SKIPF( IUNIT, 1) 

CALL  MTUAIT(IER) 

GO  TO  10 

C  READ  IN  RECORD  FROM  TAPE 

30  DO  100  1=1.30 
IMAX( I ) =0 

CALL  MTREAD( IBUF.512. IER.0. IUNIT.NUMB) 

C  FIND  MINIMUM  VALUE  FOR  EACH  PHONE 

DO  40  J-1.512 
IF ( IBUF ( J) . GE.0)GO  TO  40 
1TEST=IABS( IBUF(J) ) 

IF ( ITEST.GT. IMAX( I) ) IPEAK-J 
IF ( ITEST.GT. IMAX( I) ) IMAX( I ) ■ ITEST 
40  CONTINUE 

C  FIND  POS  MAX  FOR  FILE  MAXVL 

IMAXKD-0 
DO  45  J-1.512 

IF ( IBUF ( J) .GT. I MAXI ( I ) ) IMAX1 ( I) =IBUF(J) 

45  CONTINUE 

C  SET  DETECTION  LEVEL  AND  FIND  1ST  NEGATIVE  BELOW  LEVEL 

IDTCT-  -ICONS*(IMAX(I)/10) 

IF( IPEAK.LT. 101 )URITE(LUOUT, 48) I 

48  FORMATdX. 'PEAK  VALUE  FOUND  WITHIN  1ST  IPO  SAMPLES  PHONE:’.  12) 
IF( IPEAK.LT. 101) IPEAK-101 
DO  50  J  =  IPEAK- 100. IPEAK 
!F(  IBUF ( J) .LT. IDTCT) GO  TO  60 
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50  CONTINUE 

IFCJ.EQ.512)URITE( LUOUT,  55 ) I 

55  FORMAT! IX,  'NO  NEGATIVE  VALUE  BELOU  THRESHHOLD  -  PHONE  *.I2) 
GO  TO  100 

C  FIND  1ST  NEGATIVE  PEAK 

60  DO  70  K  =  1  -  5  12 

IFCJ+K.GT.5123GO  TO  80 
IF C IBUF C J+K-l) .LT. IBUFC J+K) )GO  TO  80 
70  CONTINUE 
80  IPOSC I) -J+K-l 
C  FIND  ZERO  CROSSING 

L-IPOSCI) 

90  L=L+1 

IF (L.GT.500)GO  TO  100 
IF C IBUF CL) .LT.0)GO  TO  90 
IPOS1 C  I)  =L 
IZERP( I ) = IBUF CL) 

IZERNC I ) =  IBUFCL-1 ) 

100  CONTINUE 

CHECK  SEQUENCE  FOR  BAD  MAX  OR  MIN 
IF  FIND  ONE  GO  BACK  OR  FORWARD  1I*N  SAMPLES 
DO  140  1=1.27 
ICHK( I) =0 
140  CONTINUE 

DO  150  1=2.26 
ITRYl-IPOSCI-l)-IPOS(I) 

ITRY2- IPOSC I)- IPOS ( 1+1) 

IF ( ITRV1 .GE.0. AND. ITRY2.GE.0)GO  TO  150 
IF ( ITRY1 .LE .0 . AND . ITRY2.LE.0)GO  TO  150 
IF ( IABS ( ITRY1) .LE. 10. OR. IABS( ITRY2) .LE. 10) GO  TO  150 

iPos<i)=(iPosn-n+iPosn+i))/2 

ICHK C I )  =  1 
150  CONTINUE 
C  SEE  IF  ANY  LOCAL  MINIMA  CHANGED 

DO  152  1=1,27 
IFCICHKU)  .NE.0)GO  TO  154 
152  CONTINUE 
GO  TO  165 

C  REPOSITION  TAPE  TO  HEADER  RECORD 

154  CALL  SKIPF ( IUNIT.0) 

CALL  MTUAITCIER) 

C  REFIND  LOCAL  MINIMA  FOR  THOSE  CHANGED  ABOVE 

DO  160  1-1.27 
CALL  SKIPBC IUNIT. 1) 

CALL  MTtJAIT(IER) 

IFCICHKCI)  .EQ.0)GO  TO  160 
JTEMP-1 

TYPE  REFIND  MINIMA  PHONE  *  '.I 
CALL  MTREADC IBUF. 512. IER.0. IUNIT. NUMB) 

155  IX- IBUF ( IPOSC I) ) 

IXM1-IBUF ( IPOSC I)-l) 

IXP1 » IBUF ( IPOSC I )+l ) 

IFCIX.LE.  IXM1.AND.  IX. LE.  IXPDGO  TO  159 
IF ( IX. GT. IXM1) IPOSC I) -IPOSC I) -1 
IFCIX.LT. IXM1.AND. IX.GT. IXP1) IPOSC I) - IPOSC I)+l 
IFC  IPOSC  I) .LE.50.OR. IPOSC I) .GE.450)GO  TO  159 
TYPE  PHONE  *  '.I,'  ITERATION  *  ' , JTEMP 
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JTEMP=JTEMP+1 
IFUTEMP.GT.  100)GO  TO  159 
GO  TO  155 

C  REFIND  ZERO  CROSSING 

159  L-IPOSU) 

158  L=L+1 

IF (L.GT.500)GO  TO  160 
IF(IBUF(L) .LT.0)GO  TO  158 
IPOS  1 ( I) =L 
IZERP !  I)  = IBUF !L) 

IZERN( I) =I8UF (L-l) 

160  CONTINUE 

C  CONVERT  SEQUENCE  FROM  POSITIONS  TO  ARRIVAL  TIMES 

165  1SEC1=ITIM!4)-15*! ITIM!4)/15) 

SEC=FLOAT ( ISEC1 )+FLOAT! ITIM!5> )/FLOATC ITIM(6) ) 

DO  170  1=1.29 

TDEL 1 “FLOAT ( IPOS  1 C I) - 1) /409G . 

IF ( IZERN! I) . EQ.0.AND. IZERP! I) . EQ.0) GO  TO  170 
TDEL2=FL0AT!-IZERN! I) )/!FLOAT! IZERP! I) ) -FLOAT! IZERN(I)) ) 
TDEL2-TDEL2/4096. 

TDEL ( I ) =SEC+TDEL 1 +TDEL2 
170  CONTINUE 
C  DEFINE  ITIM2 

DO  175  J=l,6 
ITIM2!J)=ITIM!J) 

175  CONTINUE 

C  WRITE  ONE  OUTPUT  RECORD  TO  LP  EVERY  TWO  MINUTES 

ITST=ITIM!3)-2*! ITIM!3)^2) 

IF( ITST.NE.0.OR. ITIMI4) .GE. 15)G0  TO  985 
WRITE (2. 994) !ITIM!J),J=1»4) 

994  FORMAT( IX, 'START  TIME:  '.414) 

WRITEC2.991) 

99 1  FORMAT ( 1 X, 2 ( ' PHONE ' . 4X, ' MAX' , 2X, ' NEGPK ' . 2X, ' ZERO+' , 2X, ' +  VAL ' , 

1  2X, ' -  VAL'.2X,'TM  DELAY' ,6X) ) 

DO  996  1-1.15 

WRITE (2. 995) I, IMAX( I) , IPOSC I ) . IPOS  I ( I) . IZERP( I) . IZERNC I) . TDEL ( I) 
1. 1  +  15.  IMAXU  +  15).  IPOS (1  +  15).  IPOSK  1+15).  IZERP (1+15).  IZERNC  1+15), 

2  TDELCI+15) 

995  FORMAT ! 1X.2C6C I5.2X) .F8.6.6X) ) 

996  CONTINUE 
WRITE(2.990) 

990  FORMAT  C IX. '  ') 

985  WRITE! 1' IPNT1) (IARRl(J). J-1.S0) 

WRITEC3' IPNT0) ( IARR0(J) . J-1.60) 

C  SKIP  OVER  EOF 

CALL  SKIPFdUNIT.  1) 

CALL  MTWAIT(IER) 

C  GO  TO  START  OF  NEXT  RECORD 

GO  TO  10 

C  TYPE  READ  ERROR  MESSAGE  ON  TERMINAL 

110  TYPE  *.'  READ  ERROR  * 

C  CLOSE  FILES  AND  SPOOL  OUTPUT 

99  CALL  CLOSE! 1) 

CALL  CLOSE! 3) 

CALL  SPOOL! 2) 

END 
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C  PROGRAM  MRGTM 
C 

C  THIS  FILE  WILL  EDIT  THE  THREE  FILES  TIMDL.TPN 

C  N“l,2,3  FOR  BAD  POINTS  AND  WILL  MERGE  THE  RESULT 

C  INTO  FILE  TIMDL. ALL.  FORMAT  FOR  TIMDL. ALL: 

C  UORD  #  1  RECORD  * 

C  2-7  TIME  (DAY.HR.MIN.SEC.NUMER.DENOM) 

C  8-61  TIME  DELAYS 

C 

C  PROGRAMMER:  L.  ROSENBLUM  VERSION:  APRIL  23.  1980 

C 

DIMENSION  IBUF (60). IBUF 1(61) 

EQUIVALENCE  ( IBUF 1(1). IREC) . ( IBUF 1(2), IBUF( 1)  ) 

IREC-1 

CALL  ASSIGNC2.  'SY: TIMDL. ALL' ) 

DEFINE  FILE  2  ( 1289.61.U.  IPNT1) 

IPNT1-1 

DO  100  ITP-1,3 

IF C ITP.EQ. 1 ) CALL  ASSIGN( 1. ' SY: TIMDL. TP  1' ) 

IF( ITP.EQ.23CALL  ASSIGN(1.'SY:TIMDL.TP2') 

IF ( ITP . EQ . 3) CALL  ASS IGN ( 1 . ' SY: TIMDL . TP3' ) 

DEFINE  FILE  1 (500, 60. U, IPNTR) 

IPNTR*1 

IF (ITP.EQ. 1)NUMB*396 
IF ( ITP .EQ .2) NUMB =444 
IF ( ITP.EQ. 3) NUMB =449 
DO  50  1*1, NUMB 

READ(1'  IPNTR, ERR-99)  ( IBUF(J) , J-1,60) 

IF (IBUF (5) . LT.22.0R. IBUFC5) .GT. 27) TYPE  *. (IBUFCJ) , J-1.6) 

IF ( IBUF (53 .LT.22.0R. IBUF (5) . GT.27)GO  TO  50 

IPNT1»IREC 

UR  I TE  ( 2 ' IPNT1) ( IBUF1 (J),J*1,61) 

IREC* IREC+1 
50  CONTINUE 
99  CALL  CLOSE ( 1) 

100  CONTINUE 

CALL  CLOSE (2) 

END 


C  PROGRAM  ARSH7 

C 

C  THIS  PROGRAM  DILL  COMPUTE  X.Y  COORDINATES  FOR  THE  ARRAY 

C  POSITIONS  AT  GIVEN  TIMES  AND  STORE  THE  RESULTING 

C  POSITIONS  IN  THE  FILE  XYPOS.ALL.  TIME  DELAYS  ARE  READ 

C  FROM  THE  FILE  TIMDL. ALL.  THIS  PROGRAM  IS  A  MODIFICATION 

C  OF  THE  UATERFALL  PLOT  PROGRAM  ARSH9. 

C 

C  FILE  TIMDL. ALL  HAS  THE  FOLLOWING  RECORD  STRUCTURE 

C  UORD  1  RECORD  NUMBER 

C  WORDS  2-7  TIME(DAY.HR.MIN.SEC.NUMER.DENOM) 

C  WORDS  8-61  27  TIME  DELAYS 

C 

C  FILE  XYPOS.ALL  HAS  THE  FOLLOWING  RECORD  FORMAT 

C  WORD  1  RECORD  * 

C  WORDS  2-7  T I ME ( DAY, HR, MIN. SEC, NUMER , DENOM) 

C  WORDS  8-61  27  X-POSITIONS 

C  WORDS  62-115  27  Y-POSITIONS 
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c 

C  PROGRAMMER:  L. ROSENBLUM  VERSION:  5-5-80 

C 

IMPLICIT  DOUBLE  PRECISION  (A-H.O-Z) 

REAL  TDEL.X2. Y2 

DIMENSION  D (27) , A(27) .THETA (26) ,THET1 (26) , PHI (2?) 

DIMENSION  PS  I (27), PH  1 1(2?) ,X(27) , Y(27) ,X1 (2?) , Y1 (2?) 

DIMENS  ION  ITIM(6) , TDEL (27) , IBUF (6 1 ) .PHNDL (30) , 1BUF 1 ( 1 15) , X2(27) 
DIMENSION  Y2C27) 

EQUIVALENCE  ( IBUF (2) . ITIM( 1) ) . ( IBUF(8) , TDEL( 1 ) ) . ( IBUF 1 (8) ,X2( 1) ) - 
1  ( IBUF 1 (62) , Y2( 1) ) 

DATA  SNDSP.DBET,A(1).PHI(1).PHI(2)/1522.,20..0. ,0.,0./ 

DATA  PI, YMAX1.SYSDL.PLTY/3. 1415926.5.0,-17.0,  1.0/ 

DATA  PHNDL/.080.- .007,- .0 15. .034,-. 146. .088. -.031, .020. -.002. 

1  -. 107, -.080. -.042. -.012. . 141,. 011. .094. .039,-. 108, 

2  .090, -.032, . 123.-. 106. .036. -.033. . 140,-. 133, -.029, 

3  .0.0. 0,0.0/ 

A (2) -DBET 

TYPE*. '  ENTER  START  RECORD  *' 

ACCEPT  *. INUM0 

TYPE  *. '  ENTER  #  OF  RECORDS  FOR  PROCESSING  ' 

ACCEPT  *. INUM1 
C  SET  UP  FILES 

CALL  ASSIGN( 10, 'SY:TIMDL. ALL' ) 

DEFINE  FILE  10  (  1289, 61. U.  IPNTT?) 

CALL  ASSIGN< 1 1. 'SY:XYPOS. ALL' ) 

DEFINE  FILE  11  ( 1275. 1 15, U, IPNT1) 

IPNTR= INUM0 
IPNT1-1 

DO  500  II-1.INUM1 

r  rpat^  a  RFrnpn 

READ( 10' IPNTR.ERR-99) ( IBUF ( J) . J»1 .61) 

C  DO  NOT  PROCESS  BAD  RECORDS 

IF ( 1 1 . EQ . 68 . OR .  1 1 . EQ . 178. OR . 1 1 .EQ . 191 .OR. 1 1 .EQ.676.0R. 

1  II. EG. 831. OR. II. EQ. 1206)GO  TO  500 

C  ADJUST  FOR  SYSTEM  DELAYS  AND  SINGLE  PHONE  DELAYS 

ITST1-I I-(50*( 1 1/50) ) 

IF ( ITST1 .EQ .0) TYPE  *. '  REC  *  '.II 
DO  27  1-1.27 

C  ADD  PHONE  TO  PHONE  CORRECTIONS  ' 

TDELC I ) -TDEL( I) +CSYSDL/1000 .0)+(PHNDL( D/1000.0) 

27  CONTINUE 

C  CONVERT  FROM  TIME  DELAYS  TO  DISTANCES 

DO  25  1-1.27 
D ( I ) - 1522 . *TDEL ( I ) 

25  CONTINUE 
C  GET  THETA'S 

DO  20  1-1.26 

TEMP-D( I)*D( I)+D( I+1)#D( I+1)-DBET*DBET 
TEMP-TEMP/(2.0*D( I)*D( 1+1) ) 

THETA( I) -DACOS(TEMP) 

20  CONTINUE 

C  GET  SUM  OF  THETAS 

THET 1(1) -THETA ( 1 ) 

DO  30  1-2.26 

THET1 ( I ) “THETA ( I) +THETI ( I- 1) 

30  CONTINUE 
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C  GET  A'S 

DO  40  1=3,27 

A C I )  =D ( 1 ) *D ( 1 )  +D ( I) *D (  H -2 . 0*D ( 1 ) *D ( I) *C0S  CTHET1  ( I- 1) ) 

A( I) =SQRT(A( I) ) 

40  CONTINUE 
C  GET  PHI 

DO  50  1=3,27 

TEMP-A (I - 1 ) *A C I- H +A ( I ) *A C I ) -DBET*DBET 
TEMP=TEMP/C2.0*A(I-1)*A(I» 

PHI  ( I) =DACOS(TEMP3 
50  CONTINUE 
C  GET  SIGN  PHI 

DO  60  1=3,27 

PSI  ( I)  =DASIN(  (D(  I) *SIN(THET1  (I-1)))/A(I)) 

TEMP =S IN (PI -PS  I ( I) -THET1 ( 1-2) ) 
VTST=((D(n*SIN(PSI(nn/TE^F)-D(I-l) 

IFCYTST.LT, 0, 0)PHI(I)  -PHKI) 

60  CONTINUE 

C  GET  SUM  OF  PHI'S 

PHI  1 ( 1) =0.0 
DO  90  1=2,27 
PH  1 1  ( I )  =PH  1 1  ( I  - ! }  +PH  Id) 

90  CONTINUE 

C  PRINT  OUTPUT 

DO  70  1=1,27 
OUTl  =  (180.0*PHI(n)/'PI 
0UT2  =  ( 180 . 0*PH 1 1  ( I) ) /P I 
70  CONTINUE 

DO  105  1=1.27 
X ( I ) =A ( I ) *COS ( PH  1 1 ( I ) ) 

Y(I)=A(I)*SIN(PHI1(I)) 

105  CONTINUE 

C  DO  LEAST  SQUARES  FIT  FOR  Y  AND  GET  SLOPE 

XAVG-0.0 
YAVG=0 . 0 
XIXI=0.0 
XI YI =0 . 0 
DO  150  1=1,27 
XAVG=XAVG+X( I ) 

YAVG=YAVG+Y( I) 

XIXI =XIXI+X( I ) I ) 

XI YI =XIYI+X( I) *Y( I ) 

150  CONTINUE 

XAVG=XAVG/27.0 

YAVG=YAVG/27.0 

XIXI=XIXI/27.0 

XIYI-XIYI/27.0 

SLOPE* (XI YI-XAVG*YAVG)/(XIXI-XAVG*XAVG) 

ANGLE =DAS lN(SLOPE/(SLOPE*SLOPE+l ) ) 

YINT=YAVG-SLOPE*XAVG 

C  ROTATE  THRU  ANGLE  'ANGLE';  X1(I),Y1(I)  =  ROTATED  PTS. 

C  TRANSLATE  (0.YINT)  TO  ORIGIN. 

DO  180  1=1.27 

X1(I)=X( I ) *COS ( ANGLE) +(Y(I)-Y1NT)*S IN (ANGLE) 

YI ( I) =-X( I ) *S IN (ANGLE) +(Y( I ) -YINT) *COS ( ANGLE) 

180  CONTINUE 
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C  STORE  AWAY  XY  COORDS  IN  XYPOS.ALL 

DO  185  1=1,7 
IBUF1C  I) =IBUF(I) 

185  CONTINUE 

DO  188  1=1.27 
X2( I) =X1 ( I) 

Y2(1)=Y1(I) 

188  CONTINUE 

IJRITE  til’  IPNT 1 .  ERR  =99)  ( IBUF  1(J),J=1,1 151 
500  CONTINUE 
C  CLOSE  FILES 

99  CALL  CLOSE (11) 

CALL  CLOSE (10) 

END 


PROGRAM  XYPL1 

THIS  PROGRAM  WILL  PRODUCE  WATERFALL  PLOTS  USING' THE  XY  POSITIONS 
PRODUCED  BY  ARSH7  AND  STORED  IN  THE  FILE  XYPOS.ALL. 

THIS  VERSION  OF  XYPLT  (JILL  PRODUCE  AVERAGED  6JATERFALL  PLOTS. 

PROGRAMMER:  L.  ROSENBLUM  VERSION  g/4/80 

DIMENSION  IBUF ( 115). ITIMC6) .XI (27) . Y1 (27) . YMEANC27) 

EOU I VALENCE  ( IBUF (2) . ITIMf 1 ) ) . f IBUF (8) , XI ( 1 ) ) . ( IBUF (62) . Y1 (1) ) 

DATA  YMAX1.PLTY/5. 0,1.0/ 

C  SFT  UP  FILE 

CALL  ASSIGN! 10, ' SY: XYPOS . ALL' ) 

DEFINE  FILE  10< 1275.  1  15, U. IPNTR) 

C  GET  LIMITS  FOR  DATA 

TYPE  *, ’  ENTER  START  RECORD  NUMBER' 

ACCEPT  *. INUM0 

TYPE  *, '  ENTER  *  OF  RECORDS  TO  BE  AVERAGED' 

ACCEPT  *, INUM3 

TYPE  ENTER  NUMBER  OF  AVERAGES  ON  ONE  PLOT' 

ACCEPT  *, INUM1 

TYPE  *, '  ENTER  *  OF  WATERFALL  PLOTS' 

ACCEPT  *, INUM2 
IPNTR= INUM0 

C  DO  NUM2  WATERFALL  PLOTS 

DO  999  JJJ= 1 . INUM2 

C  DO  INUM1  AVERAGES  ON  ONE  PLOT 

DO  500  1 1-1. INUM1 
C  INITIALIZE  MEAN 

DO  200  [=1.27 
YMEAN( I) =0. 0 
200  CONTINUE 

C  GET  1NUM3  RECORDS  FOR  AVERAGE 

DO  250  KK= l , INUM3 
C  READ  A  RECORD 

READ  (10' IPNTR. ERR=99) ( IBUF(J) . J-l. 115) 

C  ADD  NEW  VALUES  TO  MEAN 

DO  210  1=1,27 
YMEAN ( I ) =YMEAN( I ) +Y1 ( I) 

210  CONTINUE 
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250  CONTINUE 

DO  260  1=1.2? 

YMEANC I) =YMEAN( I)/FL0AT( INUM3) 

260  CONTINUE 
C  SCALE  DATA 

DO  130  1=1,2? 

XI  (IJ=X1  C  I)*6. 0/600.0 
YMEAN(I)=YMEAN(I)*0.5 
130  CONTINUE 

C  INITIALIZE  AND  DRAW  AXES 

C  SKIP  IN  IT  IF  HOT  1ST  TIME  THRU 

IF(PLTY.GT.9.0)GO  TO  99 
IF(PLTY.NE. 1.0)GO  TO  138 

CALL  AXISC0. 5.  1. 0, 0. ,600.,  100. .0.. 'METERS'. 6.1) 

CALL  SYMBOL (2. 0,0. 2,0. l.'TIME:  ',0.0,6) 

CALL  SYMBOL (6. 5, 4. 0,0. 08, ' 1  INCH  =  2  METERS  '.90.0,18) 
DO  135  1=2.5 
RNUM=FLOAT ( IBUF ( I ) ) 

IF(I.EQ.5)RNUM=RNUM-1 
XI 1 =2.0+0 . 4*(FL0AT( I )  —  1 . 0) 

CALL  NUM8ERIX1 1.0. 2,0. 1.RNUM.0. 0,-0) 

135  CONTINUE 

CALL  PLOTC0.5, 1.0.-3) 

C  DRAU  PLOT 

138  DO  140  1=1.2? 

IF (YMEAN ( I ) . GT. YMAX1 ) YMEAN ( I) =YMAX1 
I F ( YME AN ( I ) . LT . - YMAX 1 ) YMEAN ( I ) =-YMAXl 
YMEAN  C I ) = YME AN ( I ) +PLTY 
IFCI.EQ. 1)CALL  PLOT (0.0, YMEANC I ) , 3) 

CALL  PL0T(X1 (I) , YMEAN ( I). 2) 

140  CONTINUE 

PLTY=PLTY+0 . 2 
500  CONTINUE 
PLTY= 1 , 0 

TYPE  *, '  PLOT  #  '.JJJ,'  COMPLETED' 

99  CALL  PLOT (0 . , 0. ,999) 

999  CONTINUE 

CALL  ENDPLT 
CALL  CLOSEC10) 

END 


PROGRAM  ARSTA 

THIS  PROGRAM  READS  THE  FILE  XYPOS.ALL  AND  COMPUTES 
THE  MEAN  AND  STANDARD  DEVIATION.  SEE  PROGRAM 
ARSH?  FOR  THE  FORMAT  OF  XYPOS. 

PROGRAMMER:  LARRY  ROSENBLUM  VERSION  5-5-80 

DIMENSION  ISUF ( 1 15) ,X(2?) . Y(2?) . ITIMEC6),YMEAN(27), IN0C2?) 
DIMENSION  YSD (2?) 

EQUIVALENCE  ( IBUF (8) ,X( 1) ) . ( IBUFC62) ,Y( D) 

TYPE  *, 'ENTER  START  RECORD  #' 

ACCEPT  *. INUM0 

TYPE  #, 'ENTER  NUMBER  OF  RECORDS  TO  BE  PROCESSED' 
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ACCEPT  *. INUM1 

CALL  ASSIGN(11.'SY:ARSTA.0UT') 

DO  ?  I »1.27 
INO  <0-1 
7  CONTINUE 

C  SET  UP  FILES 

CALL  ASSIGN (10. 'SY:XYPOS. ALL') 

DEFINE  FILE  10  ( 1275. 1 15. U. IPNTR) 

IPNTR* INUM0 
DO  10  1-1,27 
YMEAN ( I) =0.0 
10  CONTINUE 

DO  50  II-1.INUM1 
C  READ  A  RECORD 

READ (10' IPNTR. ERR-100) ( IBUF(J) . J-l. 115) 

C  STORE  START  TIME 

IF(II.NE.l)  GO  TO  15 
v  DO  12  1-1,6 

ITIMEC  I)  -  IBUFU+1) 

12  CONTINUE 

C  ADD  NED  VALUES  TO  MEAN 

15  DO  20  1*1.27 

YMEAN ( I ) -YMEAN ( I ) +Y ( I ) 

20  CONTINUE 
50  CONTINUE 
DO  60  1-1,27 

YMEAN ( I ) -YMEAN ( I ) /FLOATC INUM1 ) 

60  CONTINUE 

DO  110  1-1.27 
YSD(  I) -0.0 
110  CONTINUE 

C  NOD  COMPUTE  THE  STANDARD  DEVIATION 

100  IPNTR- 1NUM0 

DO  150  I-1.INUM1 

READ( 10' IPNTR. ERR-99) ( IBUF( J) , J-l, 1 15) 

DO  120  J-l, 27 

YSD(J) -YSD(J)+Y(J)*Y(J) 

120  CONTINUE 
150  CONTINUE 

DO  130  1-1.27 

YSD( I ) *YSD( I )/FLOAT( INUM1 ) 

130  CONTINUE 

DO  140  1-1,27 

YSD ( I ) -SORT ( YSD ( I ) -YMEAN ( I ) *YMEAN ( I ) ) 

140  CONTINUE 

UR  I TE ( 1 1 . 65 ) INUM0. INUM1 

65  FQRMAT( IX. 'START  RECORD  14, SX. 'NUMBER  OF  RECORDS', 14) 

URITE( 11.70) (INO(J).J-l.  14) 

70  FORMATdX.  14(5X.  12.  IX)) 

URITE( 1 1.75) ( YMEAN (J). J-l.  14) 

75  FORMATdX.  14(F7. 4.  IX)) 

UR ITE (11.75) (YSD(J).J-l. 14) 

URITE( 1 1.80) ( INO(J) , J-15.27) 

80  FORMATdX.  13(5X.  12.  IX)) 

UR ITE (11.85) ( YMEAN (J).  J-15.27) 

85  FORMATdX,  13CF7. 4,  IX)) 

URITEdl. 85)  (YSD(J). J-15.27) 

UR ITE  d 1 .90) 
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DO  94  J=l,2? 

YOUT-2 . 0*YMEAN ( J ) /3 . 0 

1FCJ.GT. 1 ) YOUT 1 -Y0UT-2 . 0*YMEAN ( J- 1 ) /30 . 0 
UR  I TE ( 1 1 . 93 ) J , YOUT, YOUT 1 

93  FORMATdX, 'PHONE:  ' . I2.4X. 'CORR:  ' ,F8 .5.4X. 'MOD.  C0RR:',F8.5) 

94  CONTINUE 
90  FORMAT (5X) 

99  CALL  CLOSE (10) 

CALL  SPOOL ( 1 1 ) 

C  THIS  SECTION  PLOTS  OUT  RESULTS 

CALL  PLOTS 

CALL  AXIS(0.5,0.5,0..30.,5..0..'PHONE  '.6.1) 

CALL  AXIS (0.5, 0.5. -1.0. 1 . 0, 0 . 4, 90 . 0, 'METERS' .-6. 1) 

CALL  SYMBOL (1 .8. 6.0,0. 1. 'START  TIME:  '.0.0.12) 

DO  205  1-1.4 
RNUM-FLOAT ( ITIME ( I) ) 

IF ( I . EQ.4)RNUM=RNUM-1 ,0 
XI 1 =2 .0+0 . 4*FLOAT ( I ) 

CALL  NUMBER (XI 1.6. 0,0. l.RNUM.0.0,0) 

205  CONTINUE 

CALL  SYMBOL (2. 0.5. 8,0. 1 , ' *  OF  RECS  ',0.0,10) 

RNUM=FLOAT ( INUM1) 

CALL  NUMBER(2. 6, 5. 8,0. l.RNUM.0.0,0) 

CALL  PLOT(0.5,3.0.-3) 

CALL  PLOT (6 . 0 , 0 . 0, 2) 

C  SCALE  DATA 

DO  210  1=1,22 
XPLT=6 . 0*FLOAT ( I ) /30 . 0 
YTEM1  =2 . 5*C YMEAN ( I ) -YSD (  I) ) 

YTEM2-2 . 5*( YMEAN ( I ) +YSD C I) ) 

YTEM0-2 . 5* YMEAN ( I ) -0 . 02 
XTEM0=XPLT-0 . 0 1 

CALL  SYMBOL ( XTEM0 . YTEM0 , 0 . 04, ' 0  '.0.0.2) 

CALL  PLOT  (XPLT.YTEMl.-t-3) 

CALL  PLOT  (XPLT.YTEM2.+2) 

210  CONTINUE 

CALL  PLOT(0. . 0. .  999) 

CALL  ENDPLT 
END 


PROGRAM  ARST2 

THIS  PROGRAM  (JILL  COMPUTE  SIGMA  ACROSS  THE  ARRAY  FOR  A  FIXED 
TIME:  SIGMA-SQRT( ( 1^27) *SUMMATION( (Y( J)-YAVG)**2) ) ) . 

THE  PROGRAM  THEN  PRODUCES  PLOTS  OF  EITHER  SIGMA  VS.  TIME  FOR 
USER  INPUT  AVERAGE  INTERVALS  OR  HISTOGRAM  TYPE  PLOTS  OF  BIN 
SIZE  VS.  NUMBER  OF  OCCURANCES  OF  SIGMA  UITHIN  BIN.  UHICH  PLOT 
IS  CHOSEN  IS  USER  DETERMINED,  AS  IS  THE  BIN  SIZE  FOR  THE 
HISTOGRAM  TYPE  PLOT. 

PROGRAMMER:  LARRY  ROSENBLUM  VERSION:  JULY  1.  1980 

DIMENSION  IBUF ( 1 15) » Y(27) . ITIME (6) ,SIGMA( 1275.4) . IBIN( 100) 
EQUIVALENCE  ( IBUF ( 2 ) . IT IME C I ) ) , ( IBUF (62) . Y( 1 ) ) 
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TYPE  *. 'ENTER:  1-SIGMA  VS,  TIME  2-BIN  SIZE  VS.  *  OCCURANCES* 
ACCEPT  *, IANS 
IF(IAN$.EQ.2)G0  TO  5 

TYPE  *. 'ENTER  *  OF  MINUTES  TO  AVERAGE  OVER* 

ACCEPT  *, INUM0 
GO  TO  6 

5  TYPE  *. 'ENTER  HISTOGRAM  BIN  SIZE  <0.01  TO  0.5)* 

ACCEPT  *,BIN 

IFCBIN.LT. 0.01. OR. BIN. GT.0.5)GO  TO  5 
C  SET  UP  FILES 

6  CALL  ASSIGN( 1 1. 'SY: ARST2.0UT* ) 

CALL  ASSIGN! 10. *SY:XYPQS. ALL*) 

DEFINE  FILE  10  ( 1275. 1 15. U. IPNTR) 

IPNTR-1 

^HIF  SECTION  COMPUTES  SIGMA  FOR  EACH  SHAPE 

ICNT1-0 
READ  A  RECORD 

10  READ( 10' IPNTR, ERR-99) (IBUF(J).J-l. 115) 

GET  MEAN  -  NOTE:  MEAN  SHOULD  BE  0  SINCE  DONE  BY  LS  FIT 
YMEAN-0.0 
DO  20  1-1.27 
YMEAN-YMEAN+YU) 

20  CONTINUE 

YMEAN-YMEAN/27.0 
COMPUTE  SIGMA 
SUM-0.0 
DO  30  1-1.27 

SUM-SUM+ ( Y ( 1) - YME AN ) *  <  Y  < I ) YTEAN ) 

30  CONTINUE 

IPNT1- IPNTR- ICNT1 
SIGMA( JPNT1-1, 1)«ITIME(2) 

SIGMA( IPNT1-1.2) -ITIMEC3) 

SIGMAC IPNT1- 1 .3)*ITIME(4) 

S IGMA ( IPNT1 - 1 , 4) -SORTC ( 1 . 0/27 . 0) *SUM) 

ITST- IPNTR- ( 100*( IPNTR/100) ) 

IF(ITST.EQ.0)TYPE  *. IPNTR 
GO  TO  10 
CLOSE  FILES 
99  CALL  CLOSE (10) 

INUM- IPNTR-2-rCNTl 
IF C IANS .EQ. 2)G0  TO  155 

THIS  SECTION  DOES  PLOTS  OF  SIGMA  VS.  TIME 

GET  START  TIME 

TSTRT-SIGMAC1. 1)+S IGMA (1.2) /S0.0 
C  INITIALIZE  AND  DRAU  AXES 

CALL  PLOTS 

CALL  AXIS(0.5.0.5,0..8.,1..0..*TIME  (HRS)*. 10,1) 

CALL  AXIS(0. 5. 0.5.0. ..8, .1.90.. 'SIGMA  *,-6.1) 

CALL  SYNBOL <2.7,7 .8,0.2, 'MINUTE  AVERAGES  *.0.0.16) 
RNUM-FLOATC INUM0) 

CALL  NUMBER (0.5. 7. 8, 0.2, RHUM, 0 . 0 . -0 ) 

CALL  PLOT(0.5.0.5,-3) 
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C  SCALE  AND  PLOT 

S IGAV=0 . 0 
ICNT=0 

DO  150  I-l.INUM 

Tf1TST=S  IGMA ( 1 . 1 ) +S IGMA ( 1 . 2) /60 . 0+S IGMA ( I . 3) /3600 . 0 
IF ( (TMTST-TSTRT) .GE . (FLOAT ( 1NUM0) /SB .0) )G0  TO  110 
SIGAV=SIGAV+SIGMA( 1.4) 

ICNT= ICNT+1 
GO  TO  150 

C  HAVE  1NUH0  MINUTES  OF  DATA  IN  AVERAGE 

110  S1GAV=SIGAV/FL0AT( ICNT) 

XPLT=TSTRT+FLOAT ( INUM0) /120 . 0 

YPLT=10.0*SIGAV 

TEST=ABS(XPLT-XOLD) 

IF (TEST.LE.0.33)CALL  PLOTCXPLT, YPLT.2) 

IF C TEST. GT. 0 . 33) CALL  PLOTC  XPLT. YPLT.3) 

XOLD=XPLT 
SIGPV=SIGMA( 1.4) 

ICNT= 1 

TSTRT=SIGMA(I.  l)+SIGMA(I.2)/69.0 
150  CONTINUE 

CALL  PLOTC0..0..999) 

CALL  ENDPLT 
GO  TO  199 

THIS  SECTION  DOES  HISTOGRAM  TYPE  PLOT  OF  BIN  VS  *  QCCURANCES 

INITIALIZE  IB  IN 
155  DO  1G0  1=1.100 
IB  IN  Cl) =0 
160  CONTINUE 
C  COMPUTE  *  OF  POINTS  IN  EACH  BIN 

NBIN=INT(1. 0/BIN) 

IFCFLOATCNBIN) .NE . C 1 .0/BIN) )NBIN*NBIN+1 
DO  170  1  =  1.  INUM 
INDEX=1+INT(SIGMA(I.4)/8IN) 

IB  IN ( INDEX) • IB  IN ( INDEX) +  1 
170  CONTINUE 

C  FIND  MAX  VALUE  FOR  ALL  BINS 

IMAX=0 

DO  175  I  =  1 -HBIN 

IFCIMAX.LT.  IBINCD)  IMAX=IBINCI) 

175  CONTINUE 

C  FIND  SCALE  FOR  Y  AXIS  AND  DRAU  AXES 
CALL  PLOTS 

YMAX=FL0AT(6*( l+IMAX/6) ) 

YSTEP-YMAX/6.0 

CALL  AXISC0.5.0.5.0.. 1 . .0. 2.0.0. 'SIGMA  BINS'. 10.1) 

CALL  AXIS (0.5. 0.5.0. . YMAX. YSTEP, 90.0,'#  OF  OCCURANCES  '.-16,1) 

CALL  SYMBOL (2. 0.6. 2.0. l.'BIN  SIZE  -',0.0,10) 

CALL  NUMBER (2. 5. 6. 2.0. 1. BIN. 0.0, 2) 

CALL  PLOTC0.5.0.5.-3) 

C  DRAU  PLOT 

DO  180  I-l.NBIN 

XPLTL»5.0*FLOAT(I-l)*BIN  , 

XPLTR  -5.0*FLOAT(I)*BIN 
YPLT-6 .0*FLOAT( IB  INC  I ) ) /YMAX 
CALL  PLOTCXPLTL. YPLT .2) 
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CALL  PLOT  C  XPLTR , YPLT,  2 ) 
CALL  PLOT (XPLTR, 0,2) 

180  CONTINUE 

CALL  PLOTC0. .0. ,999) 
CALL  ENDPLT 
C  SPOOL  OUTPUT 

199  CALL  SPOOL ( 1 1) 

END 


PROGRAM  BMFRM 

THIS  PROGRAM  UILL  TAKE  THE  XY  POSITION  FILE  XYPOS.ALL 
AND  TAKE  A  FOURIER  TRANSFORM  TO  PRODUCE  BEAM  PATTERNS. 

THE  RESULT  IS  PLOTTED  OUT  IN  DB'S.  CAN  PLOT  EITHER  S INGLE 
SCAN  OR  TIME  AVERAGED  BEAMPATTERNS. 

PROGRAMMER:  LARRV  ROSENBLUM  VERSION:  6-26-80 

DIMENSION  IBUF ( 1 15) , ITIM(6) . Y(27) , Y1 (60) ,H(60) , AC2S6) ,B(256) , ISC4) 
DIMENSION  C(256).CDB(256) . IBUF 1(512). IBUF2(2) .CL IN (256) , BAVG(256) 
COMPLEX  A.B.ATEMP 

EQUIVALENCE  < IBUF (62) » Y( 1) ) * ( IBUF (2) , ITIM( 1 ) ) 

EQUIVALENCE  ( IBUF 1 ( 1 ) . CL  IN ( 1 ) ) , ( IBUF2 ( 1 ) , YHORM) 

DATA  PI, SSPD, FREQ. ARDEL/3. 1415926. 1522.0, 22?. 0, 10.0/ 

DATA  PH  I , HNPUR . XAXLN , IBEGN/0.0. 2. 66666?. 6. 0. 0/ 

TYPE  *, 'ENTER  PHI  -  ANGLE  OF  ARRIVAL  IN  DEGREES' 

ACCEPT  *,PH I 
TYPE  ENTER  FREQ.' 

ACCEPT  *,FREQ 
IF (FREQ . EQ . 0 . 0) FREQ=22? . 0 
PHI1=PHI*P 1/360.0 

B  TYPE  *. 'ENTER  AVERAGING  INTERVAL' 

TYPE  *,'0  ■=  1  SCAN,  1  *  N  MIN..  2  =  ?  HOURS  (ALL  DATA)' 

ACCEPT  *, IANS 

IF ( IANS.LT. 0.OR . IANS. GT. 2) GO  TO  8 
C  SET  UP  FILES 

CALL  ASSIGN( 1 1. 'SY: BMFRM. OUT'  ) 

CALL  ASSIGN (10.  'SY’.XYPOS .  ALL'  ) 

DEFINE  FILE  10 ( 12?5, 1 15.U, IPNTR) 

IPNTR= 1 

IF(IANS.EQ.2)G0  TO  14 

IF ( IANS.EQ. 0) TYPE  *. 'ENTER  RECORD  *' 

IF (IANS.EQ.0) ACCEPT  *. IPNTR 

IFdANS.NE.  l)GO  TO  14 

TYPE  *, 'ENTER  START  TIME  (HR.. MIN.)' 

ACCEPT  *. ISTH. I STM 
ISTRT»60*ISTH+ISTM 

TYPE  *, 'ENTER  INTERVAL  LENGTX  IN  MINUTES' 

ACCEPT  *. IMIN1 
ISTP» ISTRT+IMIN 1 
C  WRITE  OUT  TYPE  OF  PLOT 

UR  I TE ( 1 1 . 333 ) PH  I . FREQ 

333  FORMATdX, 'PHI  =  '  ,F8.2.4X,'FREQ  -  ',F8.2) 

14  IF (IANS.EQ.0) UR ITE( 11.340) IPNTR 
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340  FORMAT ( IX. ' SINGLE  SCAN  -  REC  *  ’.14) 

IF ( IANS. EQ. 2) WRITE (11. 345) 

345  FORMAT (IX, 'AVERAGE  -  ALL  DATA') 

IFCIANS.EQ. 1) UR  I  TEC  1 1.350) IMIN1, ISTH,  ISTM 
350  F0RMATC1X,  14,'  MINUTE  AVG .', 4X. ' START  TIMECHR. ,MIN. ) 
GET  LINEAR  ARRAY  PATTERN  FROM  FILE  ARLIN.DAT 
CALL  ASS IGN(2, ' SY: ARL  O'. DAT'  ) 

DEFINE  FILE  2(2. 512, U, IPNT1) 

IPNT1= 1 

READC2' IPNT1) (IBUF1 (J). J-1.512) 

READC2' IPNT1) ( IBUF2C J) . J= l .2) 

CALL  CLOSE (2) 

GET  HANNING  UINDOU  COEFS  H(J)  AND  COMPUTE  LAMBDA 
DO  5  1-1,60 

TEMP«PI*(FLOAT(I)-0.5)/30.0 
H C I ) =0 . 5*( 1 . 0-COS (TEMP) ) 

5  CONTINUE 

XLArfi-SSPD/FREQ 
READ  FILE 

10  READC10' IPNTR.ERR=99) (I8UF(J).J»1, 115) 
ITST»IPNTR-20*( IPNTR/20) 

IF( ITST.EQ.0. AND. IANS. NE.0)TYPE  *,'REC  *  ' , IPNTR 
SAVE  1ST  TIME 
IF ( IBEGN.EQ. l)GO  TO  II 
DO  1?  K-1,4 
ISCK)-ITIM(K) 

17  CONTINUE 

IF ( IANS.EQ.0)URITE( 11,355) (IS(J).J=1,4) 

11  IF( IANS.EQ. 2)  IREC1=2 
IF ( IANS.NE . l)GO  TO  12 
ITST= ITIM(2)*60+ITIM(3) 

IFCITST.LT. ISTRT)GO  TO  10 

IF ( ITST. EQ. ISTRT. AND. IBEGN.EQ. 0) IREC1= IPNTR 
IF ( IBEGN.EQ.0)URITE( 1 1 . 355) (IS(J).J*I,4) 

355  FORMATC IX, 'START  TIME:  ',414) 

IBEGN=1 

IF ( ITST.GE, ISTP3GO  TO  99 

12  DO  15  1=1,256 
A(I)»C0. 0,0.0) 

B(I)=(0. 0,0.0) 

15  CONTINUE 

PUT  Y'S  INTO  CORRECT  PLACE  IN  Yl'S 
DO  20  1=1,27 
Yl(2*(I+l))=Y(r) 

20  CONTINUE 

INTERPOLATE  BETWEEN  EVEN  Yl'S 
DO  30  1=5,55,2 
Yl(I)=(Yl(I-l)+Yl(l+l))/2.0 
30  CONTINUE 

EXTERPOLATE  FIRST  3  Yl'S 
DO  40  1=1.3 

Y1 ( I) =Y1 (4) +(FLOAT(4-I)#(Yl (4)-Yl(5))) 

40  CONTINUE 

EXTERPOLATE  LAST  FOUR  PHONES 
DO  50  1=1.4 

Y1 (56+1) =Y1 (5G)+(FL0AT( I)*(Y1 (56)-Yl (55) ) ) 

50  CONTINUE 
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C  COMPUTE  COMPLEX  FFT 

DO  120  J=1.60 
C1=2.0*PI/XLAMB 
T1=FL0AT(J)*ARDEL*SIN(PH!1) 
T2=Y1(J)*C0S(PHI1) 

TEMP=C 1*(T1+T2) 

A C J ) =CEXP ( CMPLX (0 . 0 , TEMP ) ) 
A(J)=H(J)*A(J) 

120  CONTINUE 

CALL  CFT (A.256.B) 

C  REAORDER  THE  CFT  OUTPUT  ARRAY 


DO  130  1= 

1, 128 

B(  I )  =A(  1  + 

128) 

B( 1+128) = 

Ad) 

130 

CONTINUE 

DO  140  1= 

1,256 

Ad)  =8  (I) 

140 

CONTINUE 

C  NORMALIZE  THE  ARRAY  AND  TAKE  POWER 

DO  150  1=1,256 
A TEMP  =256 . 0*A (I ) 

B ( I)  =  ( ATEMP*CMPLX(REAL ( ATEMP)  , -AIMAG(ATEMP)  )  ) 

C ( I) =HNPWR#REAL (B ( I ) ) 

150  CONTINUE 

IF( IANS.£Q.0)GO  TO  152 
C  ADD  NEWEST  FFT  INTO  AVERAGE 

DO  300  1=1.256 

IF ( IPNTR.EQ . IREC 1 )8AVG( I) =8AVG( I)+Cd ) 

IF  ( IPilTR . NE .  IREC  1)  BAVG (I ) = (C(I )  /'FLOAT!  IPNTR- IREC 1+1 ) ) 

1  +( FLOAT ( IPNTR- IREC! ) /FLOAT ( IPNTR- IREC 1+1) )#8AVG( I) 

300  CONTINUE 
GO  TO  10 

C  COMPUTE  DB'S 

152  CDBMX=-999.9 
INDEX=0 

DO  160  1=1.256 

CDB ( I )  =  10 . 0+LOG 10 (C ( I )  /'YNORM) 

IF (CDB( I) . GT.CDBMX) INDEX=I 
IF (CDB ( I ) . GT. CDBMX) CDBMX=CDB ( I ) 

160  CONTINUE 

URITEC 11.330) INDEX.  CDBMX 
99  CALL  CLOSE (10) 

IF(  IANS.EQ.0)  GO  TO  400 
C  COMPUTE  DB'S  FOR  AVERAGE 

CDBMX—999.9 
INDEX=0 

DO  320  1=1.256 

CDB  ( I )  ■  10 . 0*LOG  10  (BAVG  <  I )  /'YNORM) 

IF(CDB(I) .GT.CDBMX) INDEX-I 
IF (CDB ( I ) . GT.CDBMX) CDBMX=CDB ( I ) 

320  CONTINUE 

WR ITE ( 1 1 . 330) INDEX. CDBMX 
330  FORMAT ( IX. ' INDEX  =  I4.4X.'MAX  DB  -  *.F8.3) 

C  FIND  3DB  DOWN  POINTS 

ARDEL -ARDEL*COS (PH  I ) 

400  DO  450  1=130.256 

IF (CDB ( 129)-CDB( I) .LT.3.0)GO  TO  450 

ANG 1  =  ( 1 80 . 8/P  I ) +AS IN ( (FLOAT( I- !30)*XLAMB)/'(256. 0*ARDEL ) ) 
ANG2- ( 1 80 . 0/P  I ) *AS IN ( (FLOAT( I- 129) *XLAMB) /(256 . 0*ARDEL) ) 
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RAT1=CDB( 129)-CDB( 1-1) 

RAT2=CDB( 129)-CDB( I) 

RATIO=(3.0-RAT1)/(RAT2-RAT1) 

ANG3=ANG 1+RATI0*(ANG2-ANG 1 ) 

GO  TO  455 
450  CONTINUE 
455  DO  460  1=1,128 

IF (CDB(  129) -ODBC  129-1) . LT. 3 . 0)GO  TO  460 

ANG 1  =  (- 180 . 0/P  I) *AS  IN (  (FLOAT(  1- 1 ) *XLAMB) /(256 . 0*ARDEL ) ) 

ANG2= (- 180 . 0/P  I ) *AS IN ( (FL0AT( I ) *XLAMB) /<256 . 0*ARDEL) ) 
RAT1=CDB(129)-CDB( 130-1) 

RAT2=CD8 (129) -CDB (129-1) 

RATIO= (3 .0-RAT1 ) /(RAT2-RAT1 ) 

ANG4=ANG 1+RATI0*(ANG2-ANG 1 ) 

GO  TO  465 
460  CONTINUE 
465  ANG  =  ( ANG3-ANG4) /2 . 0 
UIDTH=2.0*ANG 

UR ITE ( 1 1 , 470) ANG3, ANG4,  WIDTH 

470  FORMAT ( IX, ' 3DB  PTS* ,F8.3.2X,F8.3, '  UIDTH  '.F8.3.'  DEG.') 

C  NOU  PLOT  OUT  BEAM  PATTERN 

CALL  PLOTS 

ANGMX= (180.0/PI) *ASIN(XLAMB/(2. 0*ARDEL) ) 

XAX I N = ANGMX/ ( XAXLN/2 . 0) 

ANGM1=15.0 
XAXI 1=  5.0 

CALL  AXIS(0.5,0.5, -ANGM 1 . ANGM1 , XAX 11,0.0,' ANGLE  (DEG. ) ' . 12, 1) 

CALL  AX IS (0.5, 0.5, -50. ,0. , 10. ,90. , ' RESPONSE  (DB)  ',-14,1) 

CALL  SYMBOLd. 0,6. 0.0.1, 'PHI  =  DEG. '  .0.0, 16) 

CALL  NUMBERd.  1, 6. 0,0.1, PHI. 0. 8,-0) 

CALL  SYMBOLd. 0.5. 7, 0. 1. 'FREQ  «  HZ. ',0.0, 16) 

CALL  NUMBERd.  1.5. 7,0.  1, FREQ. 0. 0,-0) 

CALL  SYMBOL (4. 2, 6. 0,0. 1,’XXXX  ACTUAL  ARRAY  ',0.0,20) 

CALL  SYMBOL (4. 2, 5. 7, 0. 1,'  LINEAR  ARRAY  ',0.0,20) 

CALL  PL0T(4.2,5.75,3) 

CALL  PLOT (4. 55, 5. 75, 2) 

IF ( IANS.EQ.0)CALL  SYMBOL d .5.6 .6,0. 2. 'S INGLE  SCAN  BEAMPATTERN  ',0. 
1  ,23) 

IF(  IANS.NE.2)CALL  SYMBOL (2 . 6,6. 2,0. 1. '  START  TIME:  ',0.0,12) 

XL0C=3.4 

XADD=0.35 

DO  185  1=1.4 

RNUM=FLOAT ( IS (I ) ) 

IF ( IANS.EQ. 1 .AND. I .EQ.4)GO  TO  185 
XP OS  1 «XLOC+< I - 1 ) *XADD 

IF( IANS.EQ. 1) CALL  NUMBER(XPOS1.6.2.0. l.RNUM.0.0,-0) 

185  CONTINUE 

IF ( IANS.GT.0)CALL  SYMBOLd .6,6.6, .2. 'AVERAGED  BEAMPATTERN  ',0.,20) 
IFdANS.EQ.  1)CALL  SYTB0L(2.6,6.4,  .  1. '  MINUTE  AVERAGE' ,0., 20) 

XMINl-FLOAT(IMINl) 

IFdANS.EQ.  1) CALL  NUMBER(2.  1,6.4,  .  l.XMINl.0.0,-0) 

IF ( IANS.EQ. 2) CALL  SYMBOL(2 . 7,6 . 4, . 1 , 'ENTIRE  EVENT' ,0. . 12) 

CALL  PLOT(0. 5.0. 5,-3) 

C  PLOT  LINEAR  PATTERN 
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DO  190  1=1,256 
1 1  =  IABS ( 129- I ) 

ANG  =  < 180.0/P  I) *AS INC (FLOAT! II) *XLAMB)/(256.0*ARDEL)) 
IF  C I . LT . 129) ANG=-ANG 

XPLT= (XAXLN/2 .0) +CXAXLN/2 .0)*(ANG/ANGM1 ) 

1FCABSCANG)  .GT.ANGMDGO  TO  190 
YPLT=5.0+CLINC  D/10.0 
IF  CYPLT.LT.0.0) YPLT=0.0 
IFCI.EQ. 1)CALL  PLOT (XPLT. YPLT.3) 

IF  <  I . NE  .  DCALL  PLOTCXPLT,  YPLT.2) 

190  CONTINUE 

C  PLOT  ACTUAL  ARRAY 

DO  200  1=1.256 
I I=IABS<  129-1) 

ANG= ( 180 . 0/P I ) *AS IN ( (FLOATCI I )*XLAMB) /C256 . 0*ARDEL) ) 
IFCI.LT.  129)  ANG=-ANG 

XPLT= (XAXLN/2 . 0) +CXAXLN/2 . 0)*(ANG/ANGM1 ) 

IF(ABS(ANG)  .GT.ANGMDGO  TO  200 
YPLT-5.0+CDBI  D/10.0 
IF(YPLT.LT.0.0)YPLT*0.0 
YHT=0.04 

XFACT«(6.0/7.0)*(YHT/2.0) 

YFACT=YHT/2.0 

CALL  SYMBOL (XPLT-XFACT.  YPLT- YF ACT.  YHT.'X  '.0.0.2) 

200  CONTINUE 

CALL  PLOTC0. 0,0. 0.999) 

CALL  ENDPLT 
CALL  SPOOL (11) 

END 
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